make_regs_bw <- function(bw){
  mod1 = lm_robust(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data = individual_level %>% filter(abs(forcing) < bw), clusters = mkid14) %>% tidy()
  mod2 = lm_robust(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data = individual_level %>% filter(abs(forcing) < bw), clusters = mkid14) %>% tidy()
  mod3 = lm_robust(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data = individual_level %>% filter(abs(forcing) < bw), clusters = mkid14) %>% tidy()
  mod4 = lm_robust(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data = individual_level %>% filter(abs(forcing) < bw), clusters = mkid14) %>% tidy()
  
  bind_rows(mod1, mod2, mod3, mod4) %>%
    mutate(bw = bw) %>%
    filter(term == "opium_legal") %>%
    return()
}

bws <- c(2500, 5000, 7500, 10000)

plot <- 
  lapply(bws, make_regs_bw) %>%
  bind_rows(.) %>%
  mutate(bw_label = paste0((bw/1000), "km"),
         bw_label = factor(bw_label, levels = c("2.5km", "5km", "7.5km", "10km")),
         bw_label = fct_rev(bw_label)) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = bw_label, shape = bw_label, color = bw_label)) +
    geom_point(position = position_dodge(width = 0.4)) +
    geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    xlab("Effect of Opium Concession System") +
    scale_color_grey() +
    theme_bw() +
    theme(axis.ticks.y.left = element_blank(),
          panel.grid.minor = element_blank(), 
          axis.title.y = element_blank(),
          panel.grid.major.x = element_blank(),
          legend.title = element_blank(),
          axis.line.y.left = element_blank(),
          axis.line = element_line(colour = "black"),
          panel.border = element_blank(),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA))

ggsave("./_4_outputs/figures/figure_4.tiff", plot = plot,  width = 5, height = 3, dpi = 300, compression = "lzw")

